// @HEADER
//*********************************************************************//
//  SiYuan: A numerical PDE solver                                     //
//  Copyright (2022) YUAN Xi                                           //
//  This Software is released under the BSD 2-Clause license detailed  //
//  in the file "LICENSE" in the top-level SiYuan directory            //
//*********************************************************************//
// @HEADER

#ifndef _CONVECTION_IMPL_HPP
#define _CONVECTION_IMPL_HPP

#include "Panzer_Workset.hpp"
#include "Panzer_Workset_Utilities.hpp"
#include "Panzer_IntegrationRule.hpp"

namespace SiYuan {

//**********************************************************************
//  SFlux
//**********************************************************************
template<typename EvalT, typename Traits>
SFlux<EvalT, Traits> :: SFlux( const Teuchos::ParameterList& p)
{
  std::string flux_name = p.get< std::string >("Name");
  Teuchos::RCP<PHX::DataLayout> data_layout = 
    p.get< Teuchos::RCP<PHX::DataLayout> >("Data Layout");
  flux_ = PHX::MDField<ScalarT>(flux_name, data_layout);
  this->addEvaluatedField( flux_ );

  m_value = XYZLib::variable<double>(p);

  this->setName(flux_name);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void SFlux<EvalT, Traits> :: postRegistrationSetup(
  typename Traits::SetupData  sd,
  PHX::FieldManager<Traits>&  fm)
{}

//**********************************************************************
template<typename EvalT, typename Traits>
void SFlux<EvalT, Traits> :: evaluateFields(
  typename Traits::EvalData  workset)
{
  ScalarT val = m_value.fetch()[0];
  flux_.deep_copy(val);
}


//**********************************************************************
//  VectorFlux
//**********************************************************************
template<typename EvalT, typename Traits>
VectorFlux<EvalT, Traits> :: VectorFlux( const Teuchos::ParameterList& p)
{
  std::string flux_name = p.get< std::string >("Name");
  Teuchos::RCP<PHX::DataLayout> data_layout = 
    p.get< Teuchos::RCP<PHX::DataLayout> >("Data Layout");
  flux_ = PHX::MDField<ScalarT>(flux_name, data_layout);
  this->addEvaluatedField( flux_ );

//  dim = data_layout->spatial_dimension;

//  auto pl_val = params.sublist("Data");
  m_value = XYZLib::variable<double>(p);

  this->setName(flux_name);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void VectorFlux<EvalT, Traits> :: postRegistrationSetup(
  typename Traits::SetupData  sd,
  PHX::FieldManager<Traits>&  fm)
{
  using namespace PHX;
  this->utils.setFieldData(flux_);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void VectorFlux<EvalT, Traits> :: evaluateFields(
  typename Traits::EvalData  workset)
{
  auto vals = m_value.fetch();
  for(int c=0;c<flux_.extent_int(0);c++)
    for(int p=0;p<flux_.extent_int(1);p++)
      for(int d=0;d<flux_.extent_int(2);d++)
        flux_(c,p,d) = vals[d];
}


//**********************************************************************
// Convection
//**********************************************************************
template<typename EvalT, typename Traits>
Convection<EvalT, Traits> :: Convection( const Teuchos::ParameterList& p) :
  ambient_( p.get< double >("Ambient Temperature") ) ,
  coefficient_( p.get< double >("Heat Transfer Coefficient") )
{
  auto ir = p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR");
  auto Basis = p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis");

  space_dim = Basis->dimension();
  num_cub_points = Basis->numPoints();
  num_basis = Basis->cardinality();
  num_cell = Basis->numCells();
  basis_name = Basis->name();
  
  std::string dof_name = p.get<std::string>("DOF Name");
  temperature_ = PHX::MDField<const ScalarT>(dof_name, Basis->functional);
  this->addDependentField( temperature_ );

  std::string flux_name = p.get< std::string >("Name");
  Teuchos::RCP<PHX::DataLayout> data_layout = ir->dl_scalar;
  flux_ = PHX::MDField<ScalarT>(flux_name, data_layout);
  this->addEvaluatedField( flux_ );

/*  Kokkos::DynRankView<double,PHX::Device> basisRef = 
     Kokkos::DynRankView<double,PHX::Device>("basisRef",num_basis,num_cub_points);
  
  Kokkos::DynRankView<double,PHX::Device> intrpCoords =
    Kokkos::DynRankView<double,PHX::Device>("intrpCoords",num_basis ,space_dim);

  Basis->getIntrepid2Basis()->getDofCoords(intrpCoords);
  Basis->getIntrepid2Basis()->getValues(basisRef, intrpCoords, Intrepid2::OPERATOR_VALUE);

  basis = Kokkos::DynRankView<double,PHX::Device>("basis",num_cell,num_basis,num_cub_points);
  Intrepid2::FunctionSpaceTools<PHX::exec_space>::HGRADtransformVALUE(basis,basisRef);*/

//  for(int b=0; b< num_basis; b++ )
//    for( int s=0; s < num_cub_points; s++ )
//        std::cout << b << "," << s << ", " << basis(0,b,s) << " ,  " << basisRef(b,s) << std::endl;
  
 // std::string n =  flux_.fieldTag().name();
  this->setName(flux_name);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void Convection<EvalT, Traits> :: postRegistrationSetup(
  typename Traits::SetupData  sd,
  PHX::FieldManager<Traits>&  fm)
{
  this->utils.setFieldData(temperature_,fm);
//  this->utils.setFieldData(flux_,fm);
  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void Convection<EvalT, Traits> :: evaluateFields(
  typename Traits::EvalData  workset)
{
//  Kokkos::DynRankView<ScalarT,PHX::Device> TQp = 
//     Kokkos::DynRankView<ScalarT,PHX::Device>("TQp",workset.num_cells,num_cub_points);
//  TQp.deep_copy(ScalarT(0.0));
  const Teuchos::RCP<const panzer::BasisValues2<double> > bv = this->wda(workset).bases[basis_index];
  for (std::size_t cell = 0; cell < workset.num_cells; ++cell) {
    for (std::size_t qp = 0; qp < num_cub_points; ++qp) {
      ScalarT T=0.0;
      for (std::size_t basis = 0; basis < num_basis; ++basis) {
        T += temperature_(cell,basis)*bv->basis_scalar(cell,basis,qp);
      }
      flux_(cell,qp) = coefficient_* (T - ambient_);
    //  std::cout << cell << ", " << qp << ", " << T << ", " << flux_(cell,qp) << std::endl;
    }
  }
/*  TQp_.print(std::cout);std::cout << "-----\n";
  temperature_.print(std::cout);
  Intrepid2::FunctionSpaceTools<PHX::exec_space>::evaluate(TQp_.get_static_view(),
                      temperature_.get_static_view(),  basis);*/

 /* for (int cell = 0; cell < flux_.extent_int(0); ++cell)
    for (int ip = 0; ip < num_cub_points; ++ip) {
      ScalarT T = 0.0;
      for(int f =0; f<num_basis; ++f) {
    //    T += temperature_(cell, f)* basis(f, ip);
      }
      flux_(cell,ip) = coefficient_* (T - ambient_);
    //  std::cout << cell << ", " << ip << ",  " << ambient_ << std::endl;
    //  std::cout << cell << " ," << ip << ", " << flux_(cell,ip) << std::endl;
    }*/
  //  flux_.print(std::cout);
	  
}


//**********************************************************************
//  Radiate
//**********************************************************************
template<typename EvalT, typename Traits>
Radiate<EvalT, Traits> :: Radiate( const Teuchos::ParameterList& p) :
  ambient_( p.get< double >("Ambient Temperature") ) ,
  coefficient_( p.get< double >("Radiate Coefficient") )
{
  auto ir = p.get< Teuchos::RCP<panzer::IntegrationRule> >("IR");
  auto Basis = p.get< Teuchos::RCP<panzer::BasisIRLayout> >("Basis");

  space_dim = Basis->dimension();
  num_cub_points = Basis->numPoints();
  num_basis = Basis->cardinality();
  num_cell = Basis->numCells();
  basis_name = Basis->name();
  
  std::string dof_name = p.get<std::string>("DOF Name");
  temperature_ = PHX::MDField<const ScalarT>(dof_name, Basis->functional);
  this->addDependentField( temperature_ );

  std::string flux_name = p.get< std::string >("Name");
  Teuchos::RCP<PHX::DataLayout> data_layout = ir->dl_scalar;
  flux_ = PHX::MDField<ScalarT>(flux_name, data_layout);
  this->addEvaluatedField( flux_ );

  std::string n = "Radiate: " + flux_.fieldTag().name();
  this->setName(n);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void Radiate<EvalT, Traits> :: postRegistrationSetup(
  typename Traits::SetupData  sd,
  PHX::FieldManager<Traits>&  fm)
{
//  this->utils.setFieldData(temperature_,fm);
  this->utils.setFieldData(flux_,fm);
  basis_index = panzer::getBasisIndex(basis_name, (*sd.worksets_)[0], this->wda);
}

//**********************************************************************
template<typename EvalT, typename Traits>
void Radiate<EvalT, Traits> :: evaluateFields(
  typename Traits::EvalData  workset)
{
  const Teuchos::RCP<const panzer::BasisValues2<double> > bv = this->wda(workset).bases[basis_index];
  for (std::size_t cell = 0; cell < workset.num_cells; ++cell) {
    for (std::size_t qp = 0; qp < num_cub_points; ++qp) {
      ScalarT T=0.0;
      for (std::size_t basis = 0; basis < num_basis; ++basis) {
        T += temperature_(cell,basis)*bv->basis_scalar(cell,basis,qp);
      }
      flux_(cell,qp) = coefficient_* (std::pow(T,4) - std::pow(ambient_,4));
    }
  }
	  
}

//**********************************************************************

}

#endif
